RNA-Sequenzierungs-Analyse - Aufgaben

Author

Patrick Metzger - Institut für Medizinische Bioinformatik und Systemmedizin

1 Einführung

Wir bitten euch hier, die grundlegenden Schritte einer RNA-Sequenzierungs-Analyse durchzuführen. Wir stellen euch bereits die auf das menschliche Referenzgenom ausgerichteten Reads zur Verfügung, d.h. euer Ausgangspunkt sind die “Raw Counts” für alle Gene jeder Probe. Zusammen mit den Counts stellen wir euch die Annotationstabelle zur Verfügung, die die Proben-IDs mit den tatsächlichen biologischen Bedingungen verknüpft. Der Datensatz stammt aus einem wissenschaftlichen Projekt, an dem wir hier in Freiburg arbeiten. Der Datensatz enthält Expressionsdaten eines Rektumadenokarzinoms, welches Varianten in BRAFD594G, KRASG12A und TP53R175H trägt. Mit dem RNA-Seq-Experiment soll untersucht werden, welchen Einfluss die Medikamente Sorafenib und Trametinib auf das Expressionsmuster der Tumorzellen haben. Um die Vergleiche vornehme zu können wurde ebenfalls eine Kontrollprobe mit Dimethylsulfoxid (DMSO) generiert. DMSO wurde gewählt, da die Medikamente darin gelöst sind, um sie den Zellen zuverabreichen. Wir betrachten im folgenden also die biologischen Bedingungen DMSO, Sorafenib (Sora) und Trametinib (Tram). Für alle Bedingungen liegen zwei biologische Replikate vor.

Das Ziel dieser Aufgabe ist es die signifikant veränderten Gene zwischen der Kontrolle (DMSO) und den beiden Behandlungen (Sora und Tram) zu finden. Des weiteren interssiert uns auch, ob es Unterschiede zwischen den Behandlungen gibt. Dabei hilft uns in beiden Fällen die Analyse der differentiell exprimierten Gene (DEG).

Im folgenden beschreiben wir die Aufgaben und geben euch weitere Hinweise und Erklärungen.

Wir werden in den nächsten Wochen schrittweise folgende Punkte bearbeiten:

  1. Importieren der Count Daten und erstellen einer Count Matrix
  2. Verknüpfung der Daten zu tatsächlichen biologischen Bedingungen
  3. Erstellen einer Tabelle, die alle Gene und ihre zugehörigen Gen-Identifikatoren enthält.
  4. Berechnung einer Hauptkomponentenanalyse (PCA), um einen ersten Überblick über den Datensatz zu erhalten
  5. Die DEG-Analyse
  6. Vergleich zwischen den Medikamenten
  7. Gene-Set Enrichment Analyse (GSEA)
  8. Visualisieren und Exportieren der Ergebnisse

Daraus resultieren folgende Aufgaben/ Fragen für euch. Wir haben diese auf die drei verbleibenden Wochen aufgeteilt. Weitere Details findet ihr fortlaufend in diesem Dokument.

Aufgaben/ Fragen Woche 4

Vorbereitung:

  1. Importiert die Count-Daten und erstellt daraus eine Count-Matrix
  2. Erstellt eine Tabelle, die alle gemessenen Gene und all ihre verschiedenen IDs enthält.
  3. Ordnet den Proben ihre biologische Bedingung zu
  4. Erstellt das DESeq2 Objekt und normalisiert die Daten
  5. Erstellt zwei QC Abbildungen, einmal mit den Roh-Counts und einmal mit den normalisierten Count-Werten.

Aufgaben/ Fragen Woche 5

DEG-Analyse:

  1. Berechnet eine Principal Component Analysis (PCA) und stellt diese grafisch dar.
  2. Baut auf dem DESeq2 Objekt auf und führt die “Differentially Expressed Genes” (DEG) Analyse durch.
    1. Bestimmt die DEGs für die beiden Vergleiche Sorafenib vs DMSO und Trametinib vs DMSO
    2. Wie viele DEGs konntet ihr pro Vergleich identifizieren? (Cutoff padj < 0.05)
    3. Wie verteilen sich die Anzahl der DEGs auf hoch- bzw. runter-reguliert? (log2FoldChange > bzw. < 0)
    4. Stellt die Anzahlen grafisch z.B. als Barplot dar.
  3. Baut die signifikanten DEGs als Übersichtstabellen ins HTML ein.
  4. Stellt die Ergebnisse der DEG-Analyse mit Hilfe eines Volcano Plots dar.

Aufgabe/ Fragen Woche 6

Vergleich der beiden Medikamente:

  1. Findet Gemeinsamkeiten und Unterschiede zwischen den beiden Medikamenten im Vergleich zu DMSO.
  2. Wendet einen Exakten Test nach Fischer an um aus den Genen veränderte Gensets/ Signalwege abzuleiten. Diese Art der Analyse nennt man funktionelle Analyse oder auch “Gene-Set Enrichment Analysis (GSEA)”.
  3. Stellt die Ergebnisse grafisch dar.

1.1 R Pakete

Verwendete R Paket:

Code
#if (!require("BiocManager", quietly = TRUE))
#    install.packages("BiocManager")
#BiocManager::install(version = "3.21")

library(DESeq2)
library(openxlsx)
library(pheatmap)
library(org.Hs.eg.db)
library(ggplot2)
library(FactoMineR)
library(ggrepel)
library(apeglm)
library(ggrepel)
library(viridis)
library(tidyverse)
library(kableExtra)
library(pheatmap)
library(genefilter)
library(EnhancedVolcano)
library(UpSetR)
library(msigdbr)
library(doMC) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassen
library(foreach) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassen
library(BiocParallel) # falls es hiermit Probleme gibt könnt ihr das Paket auch weglassen
# optional -> increase speed of analysis due to multicore processing
register(MulticoreParam(4)) # -> diese Zeile löschen, wenn es Probleme mit den Paketen doMC, foreach und/ oder BiocParallel gab

Hinweis: Falls diese nicht installiert sind, können sie mit dem Paketmanager von bioconductor BiocManager::install() oder install.packages() installiert werden.

Code
packages <- c(
  "doMC",
  "foreach",
  "BiocParallel",
  "UpSetR",
  "genefilter",
  "DESeq2",
  "openxlsx",
  "pheatmap",
  "org.Hs.eg.db",
  "ggplot2",
  "FactoMineR",
  "ggrepel",
  "apeglm",
  "viridis",
  "tidyverse",
  "kableExtra",
  "pheatmap",
  "EnhancedVolcano",
  "msigdbr"
)
BiocManager::install(packages)

1.2 Konfiguration

Code
## general config
mainDir <- "/Users/chris/Library/CloudStorage/GoogleDrive-christoph.kovacs@gmail.com/My Drive/BIDS Bioinformatik & Systembiologie 2025/KW4 Lernaufgabe/MIRACUM_BIDS_Bioinformatik_Systembiologie_RNA_Sequenzierung_Aufgabe"
analysisDir <- file.path(mainDir, "analysis")
degDIR <- file.path(analysisDir, "DEG")
gseaDIR <- file.path(analysisDir, "GSEA")
gageDIR <- file.path(analysisDir, "GSEA", "GAGE")
dir.create(degDIR, recursive = T)
dir.create(gageDIR, recursive = T)

Setzen des Arbeitsverzeichnisses: Die Funktion setwd() setz das Arbeitsverzeichnis und getwd() liest das aktuelle Arbeitsverzeichnis aus.

Code
setwd(mainDir)
getwd()
[1] "/Users/chris/Library/CloudStorage/GoogleDrive-christoph.kovacs@gmail.com/My Drive/BIDS Bioinformatik & Systembiologie 2025/KW4 Lernaufgabe/MIRACUM_BIDS_Bioinformatik_Systembiologie_RNA_Sequenzierung_Aufgabe"

Bei nicht UNIX basierten Betriebssystemen, z.B. Microsoft Windows, müsst ihr aufpassen, da dort der Backslash als Trenner verwendet wird. Ihr könnte aber, wenn ihr nicht sicher sein, auch die Funktion file.path(..., fsep = .Platform\$file.sep) vewenden, z.B.

Code
file.path(
  "~",
  "Documents",
  "Lehre",
  "BIDS",
  "Bioinformatik und Systembiologie",
  "2024",
  "BIDS_RNA_Seq_2024",
  "MIRACUM_BIDS_Bioinformatik_Systembiologie_RNA_Sequenzierung_Aufgabe_2024"
)

1.3 Funktionen

Wir stellen euch ein paar Funktionen zur Verfügung, die euch bei der Konvertierung der verschiedenen Gen-IDs unterstützen können. Alle Konvertierungen sind in der Hauptfunktion getGeneMat() zusammengefasst. Diese nimmt als Input die ENSEMBL Gen-IDs und konvertiert diese in die anderen und fasst alles zu einem “data.frame” zusammen.

Code
ensembl2entrez <- function(ensembl) {
  entrez <- mget(as.character(ensembl), org.Hs.egENSEMBL2EG, ifnotfound = NA)
  entrez <- lapply(entrez, function(i) return(i[1]))
  return(unlist(entrez))
}

entrez2ensembl <- function(entrez) {
  esbl <- mget(as.character(entrez), org.Hs.egENSEMBL, ifnotfound = NA)
  esbl <- lapply(esbl, function(i) return(i[1]))
  return(unlist(esbl))
}

entrez2symbol <- function(entrez) {
  symbol <- mget(as.character(entrez), org.Hs.egSYMBOL, ifnotfound = NA)
  symbol <- unlist(lapply(symbol, function(i) return(i[1])))
  return(symbol)
}

entrez2genename <- function(entrez) {
  symbol <- mget(as.character(entrez), org.Hs.egGENENAME, ifnotfound = NA)
  symbol <- unlist(lapply(symbol, function(i) return(i[1])))
  return(symbol)
}

getGeneMat <- function(ensIDs) {
  geneMat <- data.frame(ENSEMBL=ensIDs)
  geneMat$ENTREZ <- ensembl2entrez(geneMat$ENSEMBL)
  idxNA <- !is.na(geneMat$ENTREZ)
  sym <- entrez2symbol(na.omit(geneMat$ENTREZ))
  genename <- entrez2genename(na.omit(geneMat$ENTREZ))
  geneMat$Symbol <- NA
  geneMat$Symbol[idxNA] <- sym
  geneMat$Genename <- NA
  geneMat$Genename[idxNA] <- genename
  rownames(geneMat) <- geneMat$ENSEMBL
  return(geneMat)
}

2 Analyse

2.1 Vorbereitung

2.1.1 Importieren der Rohdaten aus dem Alignment und der Quantifizierung

Das ausrichten der Sequenzierungsschnipsel (Alignment der Reads) wurde mit dem Progamm STAR gemacht. STAR bietet außerdem die Möglichkeit auch gleich die Qunatifizierung der Expression vorzunemhen. Dabei wurde der Parameter --quantMode GeneCount verwendet. Hierzu wurde eine sogenannte GTF/ GFF Annotationsdatei benutzt, welche die Information beinhaltet, welches Gene zu welchen chromosomalen Koordinaten gehört. Wir betrachten hier ein un-stranded-RNA-Sequenzierungs Experiment.

Die einzelnen .tab Dateien beinhalten die Counts pro Gene. Dabei gibt die 1. Spalte den ENSEMBL Gen-Identifier an, z.B. ENSG00000223972 und die 2. Spalte die entsprechenden un-stranded RNA-Seq Counts. Diese beiden Spalten brauchen wir im Folgenden. Die ersten vier Reihen geben ein paar Zusammenfassungsstatistiken über die Count Datei und werden nicht benötigt.

Erstellt die Count Matrix und die Genreferenztabelle mit allen Genen und den zugehörigen IDs. Stellt die Count-Matrix, die betrachteten Gene und die Annotation der Proben zu den biologischen Bedingungen innerhalb des HTML Dokuments dar.

Code
library(readxl)

getwd()
[1] "/Users/chris/Library/CloudStorage/GoogleDrive-christoph.kovacs@gmail.com/My Drive/BIDS Bioinformatik & Systembiologie 2025/KW4 Lernaufgabe/MIRACUM_BIDS_Bioinformatik_Systembiologie_RNA_Sequenzierung_Aufgabe"
Code
# Annotation einlesen
annotation <- 
  read_excel("targets.xlsx") %>%
  mutate(group = as.factor(group)) %>%
  column_to_rownames("label") %>%
  mutate(label = rownames(.))
annotation
Code
read_counts <- function(file) {
  read_tsv(file, skip = 4, col_names = c("GeneID", "Count1", "Count2", "Count3"), show_col_types = FALSE) %>%
    select(GeneID, Count1)
}

# Eine Liste mit Count-Daten je Datei
count_list <- map(annotation$file, read_counts)

# Count-Matrix erstellen
count_matrix <- reduce(count_list, full_join, by = "GeneID") %>%
  setNames(c("GeneID", annotation$label)) %>%
  column_to_rownames("GeneID")
count_matrix
Code
# Matrix speichern
#write.table(
#  x = count_matrix,
#  file = file.path(analysisDir, "counts.txt"),
#  sep = "\t",
#  quote = FALSE
#)

# Gene (ENSEMBL-IDs) als Referenz extrahieren
gene_reference <- getGeneMat(rownames(count_matrix))
gene_reference

2.1.2 Vorbereiten des DESeq2 Objektes

Im nächsten Abschnitt wird das DESeq2-Modell erstellt. Damit wird die Grundlage für die spätere Analyse der differentiell exprimierten Gene (DEG) gelegt. Hierbei ist das Handbuch von DESeq2 sehr hilfreich. Wir importieren die Daten auf Grundlage der erstellten Count-Matrix, analog zu Count-Matrix-Import. Nach dem Import und dem Erstellen des DESeq2 Objektes müssen wir die Rohdaten einem Vorfilter Schritt unterziehen um sehr niedrige Counts zu entfernen. Hierbei is es hilfreich sich zu überlegen, wie viele biologische Replikate wir pro Bedingung haben. Idealerweise entfernt man alle Gene, die im vorliegenden Experiment in Summer über alle Bedingungen weniger als 5-10 Counts haben.

Wie viele Gene verlieren wir aufgrund von niedriger Expression?

Code
# Erstellen des DESeq2-Objekts
dds <- DESeqDataSetFromMatrix(
  countData = count_matrix,
  colData   = annotation,
  design    = ~ group
)

# Setze untreated group
dds$group <- relevel(dds$group, ref = 'DMSO') 

# Vorfilter
keep <- rowSums(counts(dds) >= 10) >= 2
sum(keep) # keeping 18,198 records
[1] 15805
Code
sum(dim(count_matrix)[1]-sum(keep)) # discarding 39,575 records
[1] 41968
Code
dds <- dds[keep,] 

2.1.3 Normalisierung und Differentielle Expression

Nachdem wir die schwach exprimierten Gene entfernt haben können wir die Normalisierung der Daten vornehmen. Dazu müssen wir die “size factors” und die “Dispersion” bestimmen, damit wir den tatsächlichen Signifikanztest (wir wollen den Wald- Test verwenden) anwenden können. Außerdem müssen wir die normalisierten Expressionswerte extrahieren, damit wir damit z.B. eine PCA erstellen können. Dazu verwenden wir die regularized logarithm (rlog) Transformation. Die Werte können dann mit der Funktion assay() extrahiert werden.

Code
# Schätze size factors
dds <- estimateSizeFactors(dds)
sizeFactors(dds)
 1_DMSO_1  2_DMSO_2  3_Sora_1  4_Sora_2  5_Tram_1  6_Tram_2 
0.9552529 0.9770969 0.9240881 1.1048597 1.0997930 0.9844989 
Code
# Schätze Dispersion
dds <- estimateDispersions(dds)
plotDispEsts(dds)

Code
# Führe Wald Test durch
dds <- nbinomWaldTest(dds)
plotMA(dds)

Code
# rlog-Transformation, behalte Gruppeninformation
rld <- rlog(dds, blind = FALSE)
rlog_matrix <- assay(rld)
head(rlog_matrix)
                  1_DMSO_1   2_DMSO_2   3_Sora_1   4_Sora_2   5_Tram_1
ENSG00000227232  6.2826098  6.2283855  6.3679370  6.3316307  6.3757344
ENSG00000238009  0.3414577  0.3414201  0.3415136  0.3558799  0.3510621
ENSG00000233750  2.4318674  2.4257760  2.4163919  2.4074100  2.4074929
ENSG00000237683  7.0804001  6.8627312  6.9807804  6.9395316  7.0272746
ENSG00000268903 -0.3143222 -0.3198386 -0.3197793 -0.3150454 -0.3101142
ENSG00000239906  3.1359541  3.0614112  3.0570808  3.1020419  3.0734826
                  6_Tram_2
ENSG00000227232  6.3469033
ENSG00000238009  0.3574054
ENSG00000233750  2.4412742
ENSG00000237683  7.0166147
ENSG00000268903 -0.3144766
ENSG00000239906  3.1084416

2.1.4 QC Plot vor und nach der Normalisierung

Zur besseren Interpretation und dem Verständnis der Normalisierung bietet es sich an die Count-Werte der einzelnen Gene als kombinierten Boxplot für die jeweiligen Bedingungen darzustellen; jeweils vor und nach der Normalisierung. Die Konvertierung der Count-Werte in log2(count + 1) hat sich in dieser Hinsicht bewährt.

Erstellt bitte einen Boxplot mit den log-transformierten Counts vor und einen Boxplot nach der Normalisierung.

Hinweis: Normalisierte Counts erhaltet ihr nach anwenden der DESeq() Funktion und der Extraktion der Werte mit counts(ddsObject, normalized = TRUE). Generell empfielt sich ggplot2 zum Zeichnen zu verweden, da es vielfälltige Möglichkeiten bietet und dadurch sehr ansehnliche Abbildungen erstellen kann. Um ein für ggplot2 passendes data.frame zu erstellen, können Pakete wie z.B. reshape2, dplyr, etc. hilfreich sein.

Code
## Roh-Counts extrahieren
raw_counts <- counts(dds)

# Log2 transformieren
log_raw <- log2(raw_counts + 1)

# Daten für ggplot umformen
df_raw <- as.data.frame(log_raw) %>%
  rownames_to_column("gene") %>%
  pivot_longer(-gene, names_to = "sample", values_to = "log2count") %>%
left_join(as.data.frame(colData(dds)) %>%
rownames_to_column("sample"), by = "sample")

# Boxplot erstellen
plot_counts_trans <- df_raw %>%
  ggplot(aes(x = sample, y = log2count, fill = group)) +   # why not group?
  geom_boxplot(outlier.size = 0.3) +
  theme_bw() +
  labs(title = "Counts (trans) pro Bedingung",
       x = "Bedingung", y = "Counts (trans)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_counts_trans

Code
## Normalisierte Counts extrahieren
norm_counts <- counts(dds, normalized = TRUE)

# Log2 transformieren
log_norm <- log2(norm_counts + 1)

# Daten für ggplot umformen
df_norm <- as.data.frame(log_norm) %>%
  rownames_to_column("gene") %>%
  pivot_longer(-gene, names_to = "sample", values_to = "log2normCount") %>%
left_join(as.data.frame(colData(dds)) %>%
rownames_to_column("sample"), by = "sample")

# Boxplot erstellen
plot_counts_normtrans <- df_norm %>% 
  ggplot(aes(x = sample, y = log2normCount, fill = group)) +
  geom_boxplot(outlier.size = 0.3) +
  theme_bw() +
  labs(title = "Counts (norm, trans) pro Bedingung",
       x = "Bedingung", y = "Counts (norm, trans)") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
plot_counts_normtrans

2.2 PCA - Principal Component Analysis

Bevor wir mit der tatsächlichen DEG-Analyse fortfahren berechnen wir eine PCA um einen ersten Eindruck von unseren Daten zu bekommen. Dazu verwenden wir die eben erstellten, normalisierten Expressionswerte (rlog-Werte). In der PCA sollen die Bedingungen als “individuals” betrachtet und die Gene als “variables”. Im Resultat soll jede Probe dargestellt und in seine Haupkomponenten zerlegt werden. Hierzu könnt ihr z.B. das Paket FactoMineR verwenden. Bitte erstellt einen sog. “Individuals Graph” mit den Proben und den ersten beiden Hauptkomponenten als Achsen. Wenn ihr das FactoMineR Paket verwendet findet ihr diese Infos unter pca$ind$coord. Um einen visuell ansprechenderen Abbildung zu erhalten würde ich empfehlen die Abbildung wieder mit dem Paket ggplot2 zu erstellen. Mit der Funktion ggsave() könnte ihr die mit ggplot2 erstellte Abbildung sehr einfach in ein geeignetes Format, z.B. PDF, png, etc. exportieren.

Erstellt eine PCA und stellt diese als Abbildung im HTML Dokument dar.

Code
### PCA: Principal Component Analysis (linear relationship analysis)
pca <- plotPCA(rld, intgroup = "group", returnData = TRUE)
# Variance explained
percentVar <- round(100 * attr(pca, "percentVar"))

# Plot
plot_pca <- ggplot(pca, aes(PC1, PC2, color = group, label = name)) +
  geom_point(size = 3) +
  geom_text(vjust = -1.2, size = 3) +
  xlab(paste0("PC1: ", percentVar[1], "% variance")) +
  ylab(paste0("PC2: ", percentVar[2], "% variance")) +
  coord_fixed() +
  theme_minimal()
plot_pca

Code
# Interpretation
# ==============
# PCA erklärt 87% der Gesamtvarianz in 2 Dimensionen
# DMSO bildet einen Cluster, die Samples kann als (linear) homogen angenommen werden
# Sorafenib bildet einen Cluster, die Samples können daher als (linear) homogen angenommen werden
# Trametinib bildet zwei Cluster, die Samples sind daher nicht (linear) homogen

MDS

Code
### MDS: Multi-dimensional Scaling (linear relationship analysis)
# Compute distance matrix from rlog-transformed data
d <- dist(t(assay(rld)))  # transpose because samples are columns
# Perform classical MDS (multi-dimensional scaling)
mds <- cmdscale(d, k = 2)  # k = number of dimensions (2D)
# Create a data frame with sample information
mds_df <- as.data.frame(mds)
mds_df$group <- colData(rld)$group  # adjust 'group' to your grouping column
# Plot with ggplot2
plot_mds <- ggplot(mds_df, aes(x = V1, y = V2, color = group)) +
  geom_point(size = 3) +
  labs(x = "MDS1", y = "MDS2") +
  theme_minimal()
plot_mds

Code
# Interpretation
# ==============
# MDS zeigt ähnliches Bild wie PCA
# Sorafenib bildet einen Cluster, die Samples können daher als (linear) homogen angenommen werden
# DMSO und Trametinib bilden jeweils zwei Cluster, die Samples scheinen daher nicht-lineare Beziehungen zu haben, die von der MDS nicht erfasst werden

t-SNA

Performing PCA
Read the 6 x 6 data matrix successfully!
Using no_dims = 2, perplexity = 1.000000, and theta = 0.500000
Computing input similarities...
Building tree...
Done in 0.00 seconds (sparsity = 0.611111)!
Learning embedding...
Iteration 50: error is 65.453220 (50 iterations in 0.00 seconds)
Iteration 100: error is 103.638598 (50 iterations in 0.00 seconds)
Iteration 150: error is 52.926096 (50 iterations in 0.00 seconds)
Iteration 200: error is 57.076918 (50 iterations in 0.00 seconds)
Iteration 250: error is 54.445751 (50 iterations in 0.00 seconds)
Iteration 300: error is 0.322697 (50 iterations in 0.00 seconds)
Iteration 350: error is 0.138587 (50 iterations in 0.00 seconds)
Iteration 400: error is 0.109798 (50 iterations in 0.00 seconds)
Iteration 450: error is 0.109646 (50 iterations in 0.00 seconds)
Iteration 500: error is 0.096316 (50 iterations in 0.00 seconds)
Iteration 550: error is 0.108245 (50 iterations in 0.00 seconds)
Iteration 600: error is 0.088462 (50 iterations in 0.00 seconds)
Iteration 650: error is 0.108734 (50 iterations in 0.00 seconds)
Iteration 700: error is 0.105920 (50 iterations in 0.00 seconds)
Iteration 750: error is 0.109760 (50 iterations in 0.00 seconds)
Iteration 800: error is 0.108264 (50 iterations in 0.00 seconds)
Iteration 850: error is 0.109467 (50 iterations in 0.00 seconds)
Iteration 900: error is 0.109448 (50 iterations in 0.00 seconds)
Iteration 950: error is 0.109438 (50 iterations in 0.00 seconds)
Iteration 1000: error is 0.109435 (50 iterations in 0.00 seconds)
Fitting performed in 0.00 seconds.

2.3 Differentielle Expressions Analyse

Um die DEG-Analyse durchzuführen müssen wir noch definieren an welchen tatsächlichen Vergleichen wir interessiert sind. Dazu können wir uns eine sog. Kontrast-Matrix erstellen. Folgende Vergleiche sind für uns interessant:

  • Sorafenib vs DMSO und
  • Trametinib vs DMSO,

da wir verstehen wollen, welchen Einfluss die beiden Inhibitoren auf unsere Zellen haben. Für beide Vergleiche bestimmen wir die DEGs und exportieren diese bis zu einem korrigierten (FDR) p-Wert < 0.05 in eine Tabelle. Hierzu verwenden wir die Funktion lfcShrink(). Damit die Tabellen für unsere Kollaborationspartner besser verständlich werden ist es essentiell, dass die Gene mit allen IDs, hautpsächlich aber dem Symbol, in den Ergebnissen enthalten sind.

Berechnet die DEGs für die beiden Vergleiche Sorafenib vs DMSO und Trametinib vs DMSO. Erstellt Tabellen, die die signifikanten Gene bis zu einem FDR-korrigierten p-Wert < 0.05 beinhalten. Achtet hier darauf auch das Symbol in der Ergebnistabelle zu haben. Ansonsten fällt die Zuordnung der Gene schwer. Zeigt diese im HTML Dokument. Exportiert die Ergebnisse der DEG-Analyse zusätzlich als Excel-Tabellen.

Bei der Darstellung der DEG-Ergebnisse innerhalb des HTML Dokuments müsst ihr nicht alle Spalten darstellen. Dies wird schnell sehr unübersichtlich. Ich würde sagen, dass ihr maximal drei Angaben braucht: Symbol, log2FoldChange und den korrigierten p-Wert. Bei der Darstellung innerhalb des HTMLs kann die Funktion kable() hilfreich sein.

Wie viele Gene sind pro Vergleich signifikant reguliert? Wie verteilt sich die Anzahl auf hoch- bzw. runter-regulierte Gene? Stellt diese Ergebnisse anschaulich dar.

Hinweis: Hierzu könnt ihr den log2FoldChange aus der DEG-Analyse verwenden.

2.3.1 Identifizierte DEGs als Tabellen

Code
tab_deg_sora <- kable(head(deg_sora, 10), caption = "Top DEGs: Sorafenib vs DMSO")
tab_deg_sora
Top DEGs: Sorafenib vs DMSO
Symbol log2FoldChange padj
TNS4 -4.276618 0e+00
FGF19 -3.548554 6e-07
SPRY4 -3.523185 0e+00
DUSP6 -3.409679 0e+00
PCK1 2.741600 0e+00
FOSL1 -2.659953 0e+00
PPP2R2C -2.606917 0e+00
BNIP5 2.260004 0e+00
DUSP2 -2.197986 0e+00
KIAA0040 -2.185101 0e+00
Code
tab_deg_tram <- kable(head(deg_tram, 10), caption = "Top DEGs: Trametinib vs DMSO")
tab_deg_tram
Top DEGs: Trametinib vs DMSO
Symbol log2FoldChange padj
TNS4 -2.510089 0
RBBP8NL 1.500431 0
EGR1 -1.434658 0
ASCL2 1.407300 0
FOSL1 -1.344055 0
BMF 1.338612 0
CYP1A1 1.332714 0
IER3 -1.230418 0
KIAA0040 -1.157227 0
MAB21L4 1.117756 0

2.3.1.1 Sorafenib vs. DMSO

Code
#kable(df[, c("Symbol", "log2FoldChange", "padj")], row.names = FALSE, "html", #col.names = c("Symbol", "log2FC", "adj. p-value"))  %>%
#  kable_styling() %>%
#  scroll_box(width = "100%", height = "800px")

2.3.1.2 Trametinib vs. DMSO

2.3.2 Visuelle Darstellung: Volcano Plot

Zusätzlich zur tabellarischen Darstellung der DEG-Ergebnisse kann man diese auch als sog. “Volcano Plots” darstellen. Dabei werden alle Gene als Punkte mit ihrem Signifikanzwert (y-Achse) und dem log2FoldChange (x-Achse) dargestellt. Die signifikanten Gene werden dabei farblich hervorgehoben. Mit dem R Paket EnhancedVolcano haben wir gute Erfahrungen gemacht.

Erstellt jeweils einen Volcano Plot für die beiden Vergleiche.

Code
if (!requireNamespace("EnhancedVolcano", quietly = TRUE)) {
  BiocManager::install("EnhancedVolcano")
}
library(EnhancedVolcano)

2.3.2.1 Volcano Plot: Sorafenib vs. DMSO

Code
plot_volcano_sora <- EnhancedVolcano(res_sora,
    lab = gene_reference$Symbol[match(rownames(res_sora), gene_reference$ENSEMBL)],
    x = 'log2FoldChange',
    y = 'padj',
    pCutoff = 0.05,
    FCcutoff = 1,
    pointSize = 2.0,
    labSize = 3.0,
    title = 'Sorafenib vs DMSO',
    subtitle = 'Differentiell exprimierte Gene',
    caption = 'log2FC > 1, FDR < 0.05',
    legendLabels = c('NS', 'Log2FC', 'FDR', 'FDR & Log2FC'),
    legendPosition = 'right')
plot_volcano_sora

2.3.2.2 Volcano Plot: Trametinib vs DMSO

Code
plot_volcano_tram <- EnhancedVolcano(res_tram,
    lab = gene_reference$Symbol[match(rownames(res_tram), gene_reference$ENSEMBL)],
    x = 'log2FoldChange',
    y = 'padj',
    pCutoff = 0.05,
    FCcutoff = 1,
    pointSize = 2.0,
    labSize = 3.0,
    title = 'Trametinib vs DMSO',
    subtitle = 'Differentiell exprimierte Gene',
    caption = 'log2FC > 1, FDR < 0.05',
    legendLabels = c('NS', 'Log2FC', 'FDR', 'FDR & Log2FC'),
    legendPosition = 'right')
plot_volcano_tram

2.4 Vergleich der beiden Medikamente

Um die Unterschiede und Gemeinsamkeiten zwischen den beiden Medikamenten bzw. deren jeweiligen Vergleiche zu DMSO zu finden, können wir das R Paket UpSetR verwenden. Da wir beide Medikament zu DMSO verglichen haben, ist das Vorgehen sinnvoll und wir können direkt Unterschiede und Gemeinsamkeiten identifizieren. Wir nehmen den Vergleich auf Ebene der signifikant regulierten Gene, getrennt nach hoch- bzw. runter-Regulation, vor. Dazu legen wir uns zwei Listen an, analog zu Basic Usage. Eine für die hoch-regulierten und eine für die runter-regulierten Gene. Danach kann der Befehl upset(fromList(list)) verwendet werden. Um an die einzelnen Gene in den verschiedenen “Sets” zu kommen, stellen wir euch eine Funktion upSetSets(list) zur Verfügung. Als Input Parameter übergebt ihr dieser Funktion die Liste, die auch für upset() verwendet wird.

Führt den Vergleich der beiden Medikamente durch.

Code
if (!requireNamespace("UpSetR", quietly = TRUE)) {
  install.packages("UpSetR")
}
library(UpSetR)

# function to obtain genes per set
upSetSets <- function(sets){
  list_names <- names(sets)
  attach(sets,warn.conflicts = F)
  res <- lapply(1:length(list_names),function(y){
    combinations <- combn(list_names,y)
    res<-as.list(apply(combinations,2,function(x){
      if(length(x)==1){
        p <- setdiff(get(x),unlist(sapply(setdiff(list_names,x),get)))
      }

      else if(length(x) < length(list_names)){
        p <- setdiff(Reduce(intersect,lapply(x,get)),
        Reduce(union,sapply(setdiff(list_names,x),get)))
      }

      else p <- Reduce(intersect,lapply(x,get))

      if(!identical(p,character(0))) p
      else NA
    }))

    if(y==length(list_names)) {
      res[[1]] <- unlist(res);
      res<-res[1]
    }
    names(res) <- apply(combinations,2,paste,collapse="-")
    res
  })
  result <- lapply(res, function(x) x[!is.na(x)])
  result <- unlist(result, recursive = F)
  result <- lapply(result,function(x) data.frame(ID=x))
  return(result)
  detach(sets)
}

2.4.1 Gemeinsamkeiten

2.4.1.1 Durch beide Medikamente hoch-regulierte Gene

Code
up_sora <- deg_sora |> filter(log2FoldChange > 0) |> pull(Symbol)
up_tram <- deg_tram |> filter(log2FoldChange > 0) |> pull(Symbol)

up_list <- list(
  Sorafenib = up_sora,
  Trametinib = up_tram
)

2.4.1.2 Durch beide Medikamente runter-regulierte Gene

Code
down_sora <- deg_sora |> filter(log2FoldChange < 0) |> pull(Symbol)
down_tram <- deg_tram |> filter(log2FoldChange < 0) |> pull(Symbol)

down_list <- list(
  Sorafenib = down_sora,
  Trametinib = down_tram
)
Code
## Plots
# Hochregulierte Gene
plot_upset_up <- upset(fromList(up_list), 
      order.by = "freq",
      mainbar.y.label = "Gemeinsame hochregulierte Gene",
      sets.x.label = "Gene pro Medikament")
plot_upset_up

Code
# Runterregulierte Gene
plot_upset_down <- upset(fromList(down_list), 
      order.by = "freq",
      mainbar.y.label = "Gemeinsame runterregulierte Gene",
      sets.x.label = "Gene pro Medikament")
plot_upset_down

Code
tab_upset_sora_up <- head(upSetSets(up_list)$Sorafenib)  # für hochregulierte Gene
tab_upset_tram_up <- head(upSetSets(up_list)$Trametinib)  # für hochregulierte Gene

tab_upset_sora_down <- head(upSetSets(down_list)$Sorafenib)  # für runterregulierte Gene
tab_upset_tram_down <- head(upSetSets(down_list)$Trametinib)  # für runterregulierte Gene

Jetzt kennen wir die Gene, die entweder durch beiden Medikamente oder auch nur in dem einen oder dem anderen verändert sind. Aber was machen wir jetzt damit? Wir können uns z.B. der Ressource “Molecular Signatures Database” (MSigDB) bedienen.

“MSigDB is a resource of tens of thousands of annotated gene sets for use with GSEA (gene-set enrichment analysis) software”.

Die “gene-set enrichment” Analyse ist eine Berechnungsmethode, mit der festgestellt wird, ob ein a priori definierter Satz von Genen statistisch signifikante, übereinstimmende Unterschiede zwischen zwei biologischen Zuständen (z.B. Phänotypen) zeigt. Anders ausgedrückt können wir damit bestimmen, ob die Gene innerhlab der oben bestimmten “Sets” einen signifikante “Funktion” haben und daraus schlußfolgern, dass diese “Funktion” in unserem Experiment verändert ist. Für die Analyse verwenden wir den Exakten Test nach Fisher oder auch hypergeometrischer Test. Die Funktion hyperG() für den Test stellen wir euch zur Verfügung.

Code
hyperG <- function(geneSets,DEgenes,universe, cutoff=0.1, mincount=2, parallel=T, adj.P.Val = F,
                   set.size = NULL){
  #' hyperG
  #' 
  #' @description Calculates Fisher's Exact test with the specified genes and the supplied gene-sets.
  #'
  #' @param geneSets list. Gene-Set the calculation is based on, e.g. go.bp
  #' @param DEgenes character vector. Gene IDs used for testing. Same identifiers as used for the gene-sets, e.g. ENTREZ IDs.
  #' @param universe character vector. Universe gene IDs.
  #' @param cutoff numeric. Cutoff used to identify sig. pathways. Default: 0.1.
  #' @param mincount numeric. Consider only pathways which contain at least mincount genes. Default: 2
  #' @param parallel boolean. Use parallel calculation. Default: TRUE
  #' @param adj.P.Val boolean. Use adjusted p-value for significance filtering. Is always calculated.
  #' @param set.size vector. Min and max size of allowed gene-sets. Default min:10 genes and max:500 genes.
  #'  
  #' @return the significant regualted pathways.
  #' @export
  #' @importFrom foreach, doMC
  
  require(foreach)
  require(doMC)
  if(parallel){
    registerDoMC(cores=detectCores())
    cores=detectCores()
  }else{
    cores=1
  }
  if(!is.null(set.size)){
    print('Set Size Limits')
    idx <- lapply(geneSets,function(x){length(x) <= set.size[2] & length(x) >= set.size[1]})
    geneSets <- geneSets[unlist(idx)]
  }
  l <- length(setdiff(universe,DEgenes))
  DElen <- length(DEgenes)
  results <- mclapply(1:length(geneSets), function(i){
    results <- matrix(data=NA,ncol=7,nrow = 1)
    colnames(results) <- c('Term','Count','Size','p-value','adj.P.Val','odds ratio','GeneIDs')
    geneSet <- intersect(universe, geneSets[[i]])
  e <- intersect(DEgenes,geneSet)
    a <- length(e)
    b <- DElen - a
    c <- length(geneSet) - a
    d <- l - c
    contigency.matrix <- cbind(c(a,b),c(c,d))
    res <- fisher.test(contigency.matrix,alternative = 'greater')
    results[1,'Term'] <- names(geneSets)[i]
    results[1,'Count'] <- a
    results[1,'Size'] <- length(geneSets[[i]])
    results[1,'p-value'] <- res$p.value
    results[1,'odds ratio'] <- res$estimate[[1]]
    # find genes annotated in the consensus term
    if(a > 0){
      genes <- intersect(DEgenes,geneSet)
      eid <- genes
      eid <- eid[order(eid)]
      results[1,'GeneIDs'] <- paste(eid,collapse="|")
    }
    return(results)
  }, mc.cores=cores)
    
  results <- as.data.frame(do.call(rbind, results))
  for(i in c(2, 3, 4, 5)){
    results[, i] <- as.numeric(as.character(results[, i]))
  }
  
  if(nrow(results) != 1){
    results <- results[order(results[,'p-value'],decreasing = FALSE),]
  results[,'adj.P.Val'] <- p.adjust(results[,'p-value'], 'BH')
  if(adj.P.Val){
    results <- as.data.frame(subset(results,results[,'adj.P.Val']<=cutoff))
  }else{
    results <- as.data.frame(subset(results,results[,'p-value']<=cutoff))
  }
    results <- as.data.frame(subset(results,results[,'Count']>=mincount))
  }else results <- as.data.frame(results)
  
  return(results)
}

Außer der Funktion benötigen wir noch die entsprechenden Gensets/ Signalwegen, die uns interessieren. Diese können wir direkt über ein R Paket von MSigDB (msigdbr) beziehen. Auch hierzu stellen wir euch eine Funktion zur Verfügung, die das abrufen und erstellen der passenden Gensets vereinfacht. Wir laden uns damit die “Hallmark Gene Sets” direkt mit den Gen-Symbolen der Gene als IDs. Das geladene Genset ist eine Liste, die direkt als Input für die Funktion hyperG() verwendet werden können.

Code
get_geneset_ag <- function(
  species = "Homo sapiens",
  category = NULL,
  subcollection = NULL,
  format = "entrez"
) {
  require(msigdbr)
  db_df <- msigdbr(species = species ,category = category, subcollection = subcollection)
  if(format == "entrez"){
    m_list = db_df %>% split(x = as.character(.$entrez_gene), f = .$gs_name)
    for(idx in 1:length(m_list)){
      m_list[[idx]] <- unique(m_list[[idx]] )
    }
    return(m_list)
  }
  if(format == "symbol"){
    m_list = db_df %>% split(x = .$gene_symbol, f = .$gs_name)
    for(idx in 1:length(m_list)){
      m_list[[idx]] <- unique(m_list[[idx]])
    }
    return(m_list)
  }
  if(format == "df"){
    return(db_df)
  }
}
Code
# load Hallmark gene set
hallmark <- get_geneset_ag(species = "Homo sapiens", category = "H", format = "symbol")

Um die Funktion hyperG() ausführen zu können brauchen wir jetzt noch ein sog. “universe”. Dieses beinhaltet alle Gene, die in unserem Experiment enthalten sind. Das “universe” ist ein Vektor, der alle Gene als Symbole beinhaltet. Es sollten keine NAs und/ oder Duplikate enthalten sein.

Erstellt das “universe” und führt die funktionelle Analye/ Gene-set Enrichment Analyse mit der hyperG Funktion und den Hallmark Signalwegen durch

Code
# Extrahiere die ENSEMBL-IDs aus dem dds-Objekt
ensembl_ids <- rownames(dds)

# Verknüpfe mit gene_reference, um die Symbolnamen zu erhalten
universe_df <- tibble(ENSEMBL = ensembl_ids) |>
  left_join(gene_reference[, c("ENSEMBL", "Symbol")], by = "ENSEMBL")

#  Erstelle das Universe
universe <- universe_df$Symbol |> unique() |> na.omit()

2.4.1.3 Gemeinsamkeiten: Durch beide Medikamente hoch-regulierte Hallmark Gensets

Code
top_enrichment_plot <- function(ds, title) {
  top_enrichment <- ds |> 
    slice_min(order_by = `adj.P.Val`, n = 10)

  return (ggplot(top_enrichment, aes(x = reorder(Term, -`adj.P.Val`), y = -log10(`adj.P.Val`))) +
    geom_col(fill = "steelblue") +
    coord_flip() +
    labs(
      title = title,
      x = "Pathway",
      y = "-log10(adj. p-value)"
    ) +
    theme_minimal())
}

  enrichment_up_sora <- hyperG(hallmark, up_sora, universe, cutoff = 0.05, adj.P.Val = TRUE)
  plot_enrichment_up_sora <- top_enrichment_plot(enrichment_up_sora, "Top Enriched Hallmark Pathways (Sorafenib UP)")
  plot_enrichment_up_sora

Code
## Interpretation (KI-unterstützt)
## ==============
## Die Behandlung mit Sorafenib führt zu einer signifikanten Hochregulation von Genen, die mit dem Zelltod (Apoptosis, p53), zellulärem Stress (zB Uv-Antwort), Immunantwort (Interferon Gamma) und Stoffwechselprozessen (Adipogenese, Häm- und Peroxisomenstoffwechsel) assoziiert sind. Dies unterstützt die bekannte cytotoxische und wachstumshemmende Wirkung von Sorafenib auf zellulärer Ebene.

enrichment_up_tram <- hyperG(hallmark, up_tram, universe, cutoff = 0.05, adj.P.Val = TRUE) # empty???
plot_enrichment_up_tram <- top_enrichment_plot(enrichment_up_tram, "Top Enriched Hallmark Pathways (Trametinib UP)")
plot_enrichment_up_tram

Code
## Interpretation
## ==============
## ?

2.4.1.4 Gemeinsamkeiten: Durch beide Medikamente runter-regulierte Hallmark Gensets

Code
enrichment_down_sora <- hyperG(hallmark, down_sora, universe, cutoff = 0.05, adj.P.Val = TRUE)
plot_enrichment_down_sora <- top_enrichment_plot(enrichment_down_sora, "Top Enriched Hallmark Pathways (Sorafenib DOWN)")
plot_enrichment_down_sora

Code
## Interpretation (KI-unterstützt)
## ==============
## Sorafenib führt zur signifikanten Herunterregulierung von Genen, die an Zellwachstum (MYC, mTORC1, G2M), Zellzyklusprogression, Entzündungsreaktionen (TNFα/NFκB) und Stressantworten beteiligt sind. Diese Ergebnisse sprechen für eine starke antiproliferative und antiinflammatorische Wirkung des Medikaments, ergänzt durch eine Hemmung zellulärer Überlebenssignale.

enrichment_down_tram <- hyperG(hallmark, down_tram, universe, cutoff = 0.05, adj.P.Val = TRUE)
plot_enrichment_down_tram <- top_enrichment_plot(enrichment_down_tram, "Top Enriched Hallmark Pathways (Trametinib DOWN)")
plot_enrichment_down_tram

Code
## Interpretation (KI-unterstützt)
## ==============
## Trametinib führt zu einer ausgeprägten Unterdrückung entzündungsbezogener Signalwege (TNFα/NFκB, IL-2, IL-6, STAT3), Zellüberlebenssignale sowie Komponenten der MAPK-Kaskade (KRAS). Diese Ergebnisse sind kompatibel mit dem Wirkmechanismus von Trametinib als MEK-Inhibitor und deuten auf eine kombinierte hemmende Wirkung auf Proliferation, Entzündung und Stressantwort hin.

3 Anhang

Code
save.image(file = "session.RData")
sessionInfo()
R version 4.5.0 (2025-04-11)
Platform: aarch64-apple-darwin20
Running under: macOS Sequoia 15.5

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: Europe/Vienna
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] generics_0.1.4              SparseArray_1.8.0          
 [3] DESeq2_1.48.1               lattice_0.22-7             
 [5] digest_0.6.37               magrittr_2.0.3             
 [7] EnhancedVolcano_1.26.0      evaluate_1.0.3             
 [9] grid_4.5.0                  RColorBrewer_1.1-3         
[11] fastmap_1.2.0               Matrix_1.7-3               
[13] jsonlite_2.0.0              ggrepel_0.9.6              
[15] GenomeInfoDb_1.44.0         httr_1.4.7                 
[17] UCSC.utils_1.4.0            scales_1.4.0               
[19] codetools_0.2-20            abind_1.4-8                
[21] cli_3.6.5                   rlang_1.1.6                
[23] crayon_1.5.3                XVector_0.48.0             
[25] Biobase_2.68.0              DelayedArray_0.34.1        
[27] yaml_2.3.10                 S4Arrays_1.8.1             
[29] tools_4.5.0                 parallel_4.5.0             
[31] BiocParallel_1.42.1         dplyr_1.1.4                
[33] ggplot2_3.5.2               locfit_1.5-9.12            
[35] GenomeInfoDbData_1.2.14     SummarizedExperiment_1.38.1
[37] BiocGenerics_0.54.0         vctrs_0.6.5                
[39] R6_2.6.1                    matrixStats_1.5.0          
[41] stats4_4.5.0                lifecycle_1.0.4            
[43] S4Vectors_0.46.0            htmlwidgets_1.6.4          
[45] IRanges_2.42.0              pkgconfig_2.0.3            
[47] pillar_1.10.2               gtable_0.3.6               
[49] glue_1.8.0                  Rcpp_1.0.14                
[51] xfun_0.52                   tibble_3.3.0               
[53] GenomicRanges_1.60.0        tidyselect_1.2.1           
[55] MatrixGenerics_1.20.0       knitr_1.50                 
[57] dichromat_2.0-0.1           farver_2.1.2               
[59] htmltools_0.5.8.1           rmarkdown_2.29             
[61] compiler_4.5.0